Shear viscosity for a moderately dense granular binary mixture 

Vicente Garz 

Departamento de Fisica, Universidad de Extremadura, E-06071 Badajoz, Spain 

Jose Maria Montaner 
Departamento de Electronica e Ingenieria Electromecdnica, 
Universidad de Extremadura, E-06071 Badajoz, Spain 
(Dated: February 2, 2008) 

Abstract 

The shear viscosity for a moderately dense granular binary mixture of smooth hard spheres 
undergoing uniform shear flow is determined. The basis for the analysis is the Enskog kinetic 
equation, solved first analytically by the Chapman-Enskog method up to first order in the shear 
rate for unforced systems as well as for systems driven by a Gaussian thermostat. As in the elastic 
case, practical evaluation requires a Sonine polynomial approximation. In the leading order, we 
determine the shear viscosity in terms of the control parameters of the problem: solid fraction, 
composition, mass ratio, size ratio and restitution coefficients. Both kinetic and collisional transfer 
contributions to the shear viscosity are considered. To probe the accuracy of the Chapman-Enskog 
results, the Enskog equation is then numerically solved for systems driven by a Gaussian thermostat 
by means of an extension to dense gases of the well-known Direct Simulation Monte Carlo (DSMC) 
method for dilute gases. The comparison between theory and simulation shows in general an 
excellent agreement over a wide region of the parameter space. 

PACS numbers: 05.20.Dd, 45.70.Mg, 51.10.+y, 47.50. +d 



I. INTRODUCTION 



An usual way of capturing the dissipative nature of granular media is through an ide- 
alized fluid of smooth, inelastic hard spheres. Despite the simplicity of the model, it has 
been shown to be quite useful in describing the dynamics of granular materials under rapid 
flow conditions ]} m ■ The essential difference from molecular fluids is the absence of en- 
ergy conservation, leading to both obvious and subtle modifications of the Navier-Stokes 
hydrodynamic equations. Although many efforts have been made in the past few years in 
the understanding of granular fluids, the derivation of the form of the transport coefficients 
remains a topic of interest and controversy. This problem has been addressed using the in- 
elastic Boltzmann equation or its dense fluid generalization, the Enskog equation. Assuming 
the existence of a normal solution for sufficiently long space and time scales, the Chapman- 
Enskog method [3-], conveniently adapted to inelastic collisions, has been applied to get the 
Navier-Stokes transport coefficients. For a monocomponent gas at low-density, the above 
coefficients have been explicitly determined as functions of the restitution coefficient □ □ □ 



ie_ accuracy of these 
. The analysis for 



from approximate solutions of the corresponding kinetic equations. T 
approximate results has been then confirmed by computer simulations^, 
dilute gases has been also extended to finite densities in the context of the revised Enskog 
kinetic theory (RET) 8]. This hydrodynamic theory succesfully models the density and 
temperature profiles obtained in a recent experimental study of a three-dimensional system 
of mustard seeds fluidized by vertical container vibrations p]. 

The majority of the studies on granular fluids are confined to monocomponent systems, 
where the particles are of the same mass and size. However, a real granular system is always 
characterized by some degrees of polydispersity in density and size, which often leads to 
segregation of an otherwise homogeneous mixture. Needless to say, the analysis of transport 
for multicomponent systems is much more involved than for a monocomponent gas. Not 
only the number of transport coefficients is higher but also they are functions of parameters 
such as the mole fractions, the mass ratios, the size ratios and the restitution coefficients. 
For this reason, m08t of the previous 8 tu dl es Q ate restricted to nearly elastic spheres. 
In addition, they usually assume energy equipartition so that the partial temperatures % 
are made equal to the global granular temperature T. Nevertheless, recent experiments of 
vibrated mixtures in three and two dimensions clearly show the breakdown of energy 
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equipartition. Related findings have also been reported by using kinetic theory tools |13 
and computer simulations (la. Q]. To the best of our knowledge, the only kinetic theory 
derivation of hydrodynamics for a granular binary mixture at low-density which takes into 
account nonequipartition of granular energy has been made by Garzo and Dufty^]. They 
solved the Boltzmann equation by applying the Chapman- Enskog method to obtain the 
Navier-Stokes equations and detailed expressions for the transport coefficients. In the case 
of the shear viscosity, the reliability of the kinetic theory predictions have also been assessed 
[3| in a wide parameter space by comparing those predictions with the results obtained 
from a numerical solution of the Boltzmann equation by means of the Direct Simulation 
Monte Carlo (DSMC) method[19]. The comparison shows an excellent agreement between 
theory and simulation. 

The objective here is to extend the analysis carried out in Ref. |l8( for the shear viscosity 
to higher densities by using the RET. The RET for elastic collisions (2^ is known to be 
an accurate theory over the entire fluid domain. Its generalization to inelastic collisions 
is straightforward (see, for example, Ref. 21() and the Chapman- Enskog method can be 
applied to obtain the transport coefficients. However, the derivation of the hydrodynamic 
equations for a binary mixture described by the RET is more complicated than in the case of 
the Boltzmann equation, due mainly to the technical difficulties associated with the spatial 
dependence of the pair correlation function. To simplify this analysis, here attention is 
restricted to the special hydrodynamic state of uniform shear flow (USF). At a macroscopic 
level, this state is characterized by constant partial densities n^, uniform temperature T and 
a linear flow velocity profile Uj = ay% a being the constant shear rate. For this particular 
problem the RET reduces to the original phenomonological kinetic theory proposed by 



Enskog 



221 ] . We solve the Enskog equation up to first order in the shear rate and evaluate 



both kinetic and collisional transfer contributions to the shear viscosity. This transport 
coefficient is expressed in terms of the solution of a set of coupled linear integral equations, 
which are then solved approximately (first Sonine polynomial a ppr oximation) just as in the 
case of elastic collisions. As done in the low-density analysis [la], the Sonine solution is 
compared with a numerical solution of the RET by using the Enskog Simulation Monte 
Carlo (ESMC) method H, which is an extension to the Enskog equation of the well-known 
DSMC methodQ- 

In a molecular fluid under USF, unless a thermostating force is introduced, the tempera- 
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ture grows in time due to viscous heating. As a consequence, the average collision frequency 
u(t) oc T 1//2 (t) increases with time and the reduced shear rate a* = ajv goes to zero in the 
long time limit. This fact allows one to identify in the simulation the Navier- Stokes shear 
viscosity coefficient rj for sufficiently long times. This route has been shown to be quite 



efficient to measure r] for dilute and dense gases [23, |24| . For a granular fluid, there is an 
additional energy sink term in the balance equation for the temperature competing with 
the viscous heating term. However, if the effect of the former term is exactly compensated 
by for the action of an external driving force, the viscous heating prevails and the shear 
viscosity can be again identified in the limit a* — > 0, just as in the elastic case. This was 
the procedure followed in Ref. [18] to measure rj from the simulation in the long time limit. 
It must be noted that the value of r\ calculated in this way (driven case) not necessarily 
coincides with the value of the shear viscosity obtained in the free cooling case (unforced 
case) . 

There are several motivations for this study. First, we want to assess to what extent 
the previous results obtained for the low-density regime are indicative of what happens for 
finite densities. Second, the comparison between theory and simulation allows one to check 
the degree of reliability of the approximate Sonine solution over a wide region of parameter 
space. Finally, by extending the Boltzmann analysis to higher densities, comparison with 
molecular dynamics simulations become practical. This comparison would determine the 
validity (or limitations) of the kinetic and hydrodynamic descriptions for granular flow. Such 
a test is essential to clarify the frequently made speculation that the above descriptions of 
granular flow are limited to weak dissipation. Some previous comparisons jlfil. 1^ support 
the hydrodynamic description, beyond complications due to possible instabilities. 

The plan of the paper is as follows. In Sec. [H] we review the Enskog theory and deduce 
the associated macroscopic conservation equations. The Chapman-Enskog method is applied 
in Sec. IIHI to solve the Enskog equation in the USF state through first order in the shear 
rate. An explicit expression for the shear viscosity coefficient is obtained in Sec. IIVI by 
using a lowest order expansion in Sonine polynomials. This transport coefficient is given in 
terms of the restitution coefficients, the temperature, the solid fraction, and the parameters 
characterizing the mixture (masses, sizes, concentrations). Section IS deals with the Monte 
Carlo simulation of the Enskog equation particularized to USF. The comparison between 
theory and simulation is carried out in Sec. IVI[ while a brief discussion on the relevance of 
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the results obtained is given in Sec. I VI II 



II. ENSKOG KINETIC THEORY AND CONSERVATION LAWS 

We consider a binary mixture of smooth hard spheres of masses mi and 777.2, and di- 
ameters <7i and cr 2 . The inelasticity of collisions among all pairs is characterized by three 
independent constant coefficients of normal restitution an, a 22 , and cti 2 = 021, where aij 
is the restitution coefficient for collisions between particles of species i and j. Due to the 
intrinsic dissipative character of collisions, in order to keep the system under rapid flow 
conditions it is usual to introduce an external driving force (thermostat) which does work 
to compensate for the collisional loss of energy This mechanism of energy input (different 
from those of shear flows or flows through vertical pipes) has been used for many authors in 
the past years to stud y d ifferent problems, such as non-Gaussian properties of the velocity 

], long-range correlations [28^ . collisional statistics and short- 
*J, or _ t propel . In this paper, for simplicity, we introduce 



distribution function |2fil . 
scale structure 



a deterministic force proportional to the peculiar velocity V (Gaussian thermostat). This 
thermostat has been frequently employed in nonequilibrium molecular dynamics simula- 
tions of elastic particles j3l|. Under these conditions, the Enskog kinetic equation for the 
one-particle velocity distribution function of species % is given by 

(ft + vi-V)/ i + i^-(Vi/,) = ^^[r,v 1 |/ i (t),/ i (t)] , (1) 

1 3=1 

where the constant £ is chosen to be the same for both species. Here, Vi = — u, u being 



the flow velocity. The Enskog collision operator JfAfi, fj] is [21 



x [ a ij 2 Xij( r , r - fij)fi( r , vi ;t)fj (r - a i:j , v' 2 ;t) 

-Xij( r , r + cr ij)/i( r > vi; t)fj(r + (T i:j , v 2 ; t)} , (2) 

where = a^a, with = (<Ji + aj) /2 and a is a unit vector directed along the line of 
centers from the sphere of species i to the sphere of species j upon collision (i.e. at contact). 
In addition, is the Heaviside step function, and g = Vi — v 2 . The primes on the velocities 
denote the initial values {v' 1; v 2 } that lead to {v 1; v 2 } following a binary collision: 

v'i = vi - Hji (l + a^ 1 ) (ct ■ g)a, v 2 = V2 + Hij (l + a^. 1 ) {a • g)tr, (3) 
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where //^ = m;/ (m^ +rrij). Finally, Xijl r , r + ^ijli 77 ^}] is the equilibrium pair correlation 
function of two hard spheres, one of species i and the other of species j, at contact, i.e., when 
the distance between their centers is Oij. In the original phenomenological kinetic theory 
of Enskogj^] (which is usually referred to as the standard Enskog theory), the Xij are the 
same functions of the densities {n^} as in a fluid mixture in uniform equilibrium. Here, 

n t = J dvfi(v) (4) 

is the number density of species i. On the other hand, this choice for Xij leads to some 
inconsistencies with irreversible thermodynamics. In order to resolve it, van Beijeren and 
Ernst jstj proposed an alternative generalization to the Enskog equation for mixtures, which 
is usually referred to as the revised Enskog theory (RET). In the RET, the Xij are the same 
functionals of the densities {n^} as in a fluid in nonuniform equilibrium. This fact increases 
considerably the technical difficulties involved in the derivation of the general hydrodynamic 
equations from the RET j^l. . unless the partial densities are uniform. 

The macroscopic balance equations for the particle number of each species, the total 
momentum and the total energy follow directly from Eq. by multiplying by 1, m^v, and 
imjf 2 , respectively, integrating over v, and summing over i. They are given by 

JUt + V • fail) + ^ = , (5) 

at mi 



d 

— u + u-Vu + p- 1 V-P = 0, (6) 

at 



2 



-T + u.VT--£^ +S; (V.q + P:Vu) = -(C-?)T. (7) 

1=1 

Here, ( is the cooling rate due to inelastic collisions among all species. The flow velocity u 
and the "granular" temperature T are defined by 

2 r 

P u = ^2 dvm i v fi(v) , (8) 
i=i 

- T = E/ rfv T^ (v) ' (9) 

i=l J 
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where n = n\ + n 2 is the total number density, and p = mini + m 2 n 2 is the total mass 
density. The mass flux jj for species i relative to the local flow is given by 



ii = m i J dvV /i(v). 



(10) 



The pressure tensor P and the heat flux q have both kinetic and collisional transfer con- 
tributions, i.e., P = P k + P c and q = q k + q c . The kinetic contributions are given by 



2 r 



i=i 

2 



(11) 



(12) 



while the collisional transfer contributions to the pressure tensor and the heat flux are, 
respectively, 

2 2 



3 minij 1 + oiij 
-f ^— ' lj m» + rrij 2 

i=l j=l J 



x / d\xij[r - (1 - A)<r«,r + Ao«]/<(r - (1 - A)<r«, Vi; t) fJr + A<r«, v 2 ; t), 



(13) 



2 2 



q c M) = EE 



i=i j=i 



J mi + rrij 



"j 



Ami + m» 



x / ^Axij[r - (1 - X)a- i j,T + Xa ij )f i (T - (1 - A)<r^, Vi;t)/j(r + A<r ii ,v 2 ;t). 
Jo 

(14) 

Here, G»j = /%Vi + /ijjV 2 is the center-of-mass velocity. Finally, the cooling rate £ in Eq. 
© is 

2 2 



cm) = ^EE 



a 



2 nriimj 



i=l j=l J 

(r, r + <Tij)fi(T, vi; t)/j-(r + try, v 2 ; t). 



1 — « 2 ) / dvi / dv 2 / dSQ){cr ■ g) (<x 



(15) 



The derivation of Eqs. (fT3*j) -()15 j) is given in Appendix El The collisional transfer contribu- 
tions are due to the derealization of the colliding pair and the additional density dependence 
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of the RET. They vanish in the low density limit but dominate at high densities. In the case 
of mechanically equivalent particles (mi = m 2 , an = a 22 = olvz = Qf, <J\ — <7 2 = a, Yjj = X)i 



Eqs. ([13 |) -([T5 |) reduce to those previously obtained in the monocomponent case [21]. 

The balance equations contain the mass flux, the heat flux, and the pressure tensor 
as specific averages over the distribution functions /j. The Chapman- Enskog method [3] 
provides a solution of the RET for states with small spatial variations in the form 

/i(r,vi,t) = /.[viMr.^Tfo^uCM)]. (16) 

This means that all space and time dependence of /j(r, Vi,t) occurs entirely through a 
functional dependence on the hydrodynamic fields. Such a solution is called normal and it 
is the basis for a fluid dynamics description of granular materials. Regarding the energy 
input mechanism we see that, according to the energy balance equation (J7J), the existence of 
a driving with the choice £ = ( compensates for the cooling effect due to the inelasticity of 
collisions. In that case, the macroscopic balance equations look like those of a conventional 
mixture with elastic collisions, although the transport coefficients entering in the constitutive 
equations are in general different from those of a gas of elastic particles. However, the 
evaluation of the complete transport coefficients of the RET for a multicomponent granular 
mixture is a very hard task and here we will pay attention to the shear viscosity coefficient 
only. Specifically, this coefficient will be determined in a particular simple situation (uniform 
shear flow) where the velocity field is the only inhomogeneity present in the system. In this 
case, the Xij are uniform so that the standard and revised Enskog theories are equivalent 
in this problem. Further, the simplicity of this state allows us to check our theoretical 
predictions for the shear viscosity with those obtained from a numerical solution of the 
corresponding Enskog equation. 



III. SHEAR VISCOSITY OF A DENSE GRANULAR BINARY MIXTURE 

As said above, we want to solve the Enskog equation JIJ in the specific state of the 
uniform shear flow (USF). In this state, the partial densities and the temperature T are 
uniform, while the velocity field is due to a simple shear 

U! = u 2 = u = gm/x, a = — — = constant. (17) 



S 



The temperature changes in time due to the competition between two mechanisms: on the 
one hand, viscous heating and, on the other hand, energy dissipation in collisions. Under 
these conditions, the mass and heat fluxes vanish by symmetry reasons and the (uniform) 
pressure tensor P is the only nonzero flux of the problem. The relevant balance equation is 
that for the temperature (J7J), which reduces to 

d t T+ l-aP xy = -((-t)T. (18) 

At a microscopic level, the USF is generated by Lees-Edwards boundary conditions 
which are simply periodic boundary conditions in the local Lagrangian frame V = v — a ■ r 
and R = r — a ■ rt. Here, a is the tensor with elements a n p = aS ax 5/3 y . In terms of the above 
variables, the velocity distribution functions are uniform 35j] 

/i(r,v,*) = /i(V,t), (19) 

and the Enskog equation takes the form 

dtfi - "Vy-^fi + ■ (V/,) = 4 WfM, /;(*)] • (20) 

x i=i 

In the Lagrangian frame, the Enskog collision operator Jf^ [V|/»(t), fj(t)] becomes 

V, J). J) = alxij j dV 2 J d*G(*- S )(*- S ) 

x [a^fiiVuVfjiV'z + aa t3 a y %t) - f l {V l ,t)f J {\ 2 - aa l3 a y % t)] .(21) 

Here, we have taken into account that Xij is uniform in the USF problem. Finally, the 
expressions for the collisional transfer contribution to the pressure tensor P c and the cooling 
rate £ in the Lagrangian frame are 

pc = ^E^4(i + «,') / ^ / dv 2 [ ^e(£.g)(£. g ) 2 

i=l j=l 1 3 ' 

xaafi (Vi + aa tj a y % t) fj{V 2 , t), (22) 

^^tt ^-x«4(i - 4) / ^ / ^ / <» S(9 • g)(9 ■ g)« 

i=l j = l 3 J J J 

xfi (V x + acrijdyx, t) fj(y 2 , t). (23) 
The normal solution for the USF state adopts the form 

f i {r,v,t)=f i (V,T(t)), (24) 



i.e. all the space dependence is accounted for by the flow velocity while all the time de- 
pendence appears through the temperature. The Chapman-Enskog method provides this 
normal solution as an expansion for small spatial gradients, i.e., as a power series in the 
shear rate a: 

/< = /, <W +/, (1) + -, (25) 

(k) 

where fl is of order k in a. The time derivatives of the fields, the Enskog collision operator, 
and the pressure tensor are also expanded as 

d t =dr ) +di i) +..., j*=js ) +js ) +---, (26) 



p = p(°) + p« + . . . . (27) 

The coefficients in the time derivative expansion are identified by a representation of the 
momentum flux, the cooling rate, and the external parameter force £ in the energy balance 
equation (|18|) as a similar series through their definitions as functionals of /j. Consequently, 
the action of the operator d t is 

0(°>T=-(C< o >-£< o >)T, d t (1) T = 0, (28) 



d\ k >T = -—aP^-V - ((( fc ) - £( fc ))T, k > 2. (29) 
3n y 

Upon writing these equations we have taken into account that P^ = = £^ = 0. The 
last equality follows from the fact that the cooling rate is a scalar, and contributions to ( in 
the first order in the gradients can arise only from V • u, which is zero in the USF. 
The leading term is the solution to the nonlinear equation 

df'f?' + y ■ (v/f >) = E J\f lft\ ff\ 00) 

3=1 



where 



a 



- 2 /i (0) (Vi)/i 0) (V^-ZfCVO/f (V 2 ) . (31) 



■I.I 



Dimensional analysis requires that ^ (V) must be of the form 

ff ) \v)=n l vfUV/vo) (32) 
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where 



v 



\ 



(33) 



i=i 



is a thermal velocity defined in terms of the temperature T of the mixture. According to 
(J32J), the time derivative in (f3"U|) can be represented more usefully as 

d 



d (V f (o) = _ (C ( 0) _ e (o) )raT/ (o) = 1 (C (0) _ e ( 0)) ^_ . ( v/ < 
The Enskog equation at this order can be written finally as 

ic <o »^-(v/r)=x:jf[/r,/f»]. 



(34) 



(35) 



Therefore, Eq. (JST?)) happens to be formally identical to the one obtained in the unforced 
case (i.e., with £ = 0) [13], and consequently there is an exact correspondence between the 
homogeneous cooling state and this type of driven steady state. This is one of the advantages 
of the Gaussian thermostat. Since the distribution functions // are isotropic, the zeroth 



order pressure tensor is found from Eqs. (fTTj) and as P^J = p5 a/ 3, where the pressure p 
is 



P 



i=l i=l 3=1 1 3 J J 

x y d*<d){a-z)(S -g) 2 



2tt 



^ + — ^ XI ^jXijnirijfijiil + a y )Ti. 



(36) 



8=1 i=l J=l 

Here, we have introduced the kinetic temperatures for each species defined as 



-n;T; 



(37) 



2 J 2 

As said in the Introduction, in general the partial temperatures T, differ from the (global) 
temperature T and so the total energy is not equally distributed between both species 
(breakdown of energy equipartition). 

The analysis to first order in a is worked out in Appendix [BJ The distribution /^ obeys 
the integral equation 

2 



3=1 



(38) 
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A similar equation can be obtained for f% \ by just making the changes 1 *-* 2. The specific 
form of the linear operators £j, A4i, and Ay are also given in Appendix[Bj The contributions 
/ i ( - ' 1 and ft determine the pressure tensor to first order in the shear rate. The result 
is 

P af5 = ~ T l a (Saxdpy + 5 ay 8p x ) , (39) 

where r] is the shear viscosity coefficient. This coefficient has kinetic and collisional transfer 
contributions 

rj = r] k + r] c . (40) 

The kinetic contribution i] k is given by 

2 



1=1 

while the collisional contribution r)° is 

2 2 



rt = ~ fdVV.VyfPty), (41) 



2 2 

% EE 4x^(1 + ti+^f dVi [ dV 2 ft\^ 0) (V 2 )g 



i=l j=i 

IV. SONINE POLYNOMIAL APPROXIMATION 



(42) 



For practical purposes the integral equations (|33j) and (J3~%j) for and are solved by 
using low order truncation of expansions in a series of Sonine polynomials. The polynomials 
are defined with respect to a Gaussian weight factor whose parameters are chosen such that 
the leading term in the expansion yields the exact moments of the entire distribution with 
respect to 1, m^v and \rriiV 2 . In the leading order, the distribution $j appearing in Eq. (J52J) 
is given by 



CO - 3 e 



7T 



4 v« 



1 + I f #^* 4 - 5^V* 2 + H 



(43) 



where 1/* = V/uo, 



2 

m, 



e i ='-^ y £mj\ (44) 



7i ,=i 



and 7, = Tj/T. For elastic collisions, 7j = 1, i.e., the partial temperatures Tj coincide with 
the global temperature T. In the inelastic case, 7« 7^ 1 and presents a complex dependence 
on the parameters of the problem. The coefficients Cj (which measure the deviation of $j 
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from the reference Maxwellian) are determined consistently from the Enskog equation. The 
approximation ()43|) provides detailed predictions for the cooling rate the temperature 
ratio Ti/T 2 and the cumulants q as functions of the mass ratio, size ratio, composition, 



density, and restitution coefficients ^j]. Recently, the accuracy of this approximate solution 
has been confirmed by Monte Carlo and molecular dynamics simulations over a 
wide range of values in the parameter space. 

In the case of the distributions £ , the leading Sonine approximation is 



(i) 



af ijM ^V x V y , f i>M (V) = n l (m l /2T t ) 3 / 2 ex P (-m l V 2 /2T l ). 



n{T 2 



(45) 



By using ()45|). the partial kinetic contributions rfc to the shear viscosity can be obtained from 
Eq. (J38j) by multiplying it with miV x V y and integrating over the velocity. From dimensional 
analysis rjf oc T 1//2 and so one gets the coupled set of equations 

/ 



rii-^ (0) + C (0) ) 

t 22 -K£(°) + c {0) ) 

where 



V 



rr\/n x T 2 
ri\ln 2 T 2 




n{T 2 



dVrrnVxVyCiif, 



(46) 



(47) 



dVm^VyMiU) 



(i + h 



(48) 



A 



dVm^VyK^jf] 



(49) 



The integrals ()49j) are evaluated in Appendix while the collision integrals ()47|) and 
were already evaluated in the Boltzmann limit (except for the factor Xij)- The explicit 
form of these integrals are also quoted in Appendix The solution of ()46|) with the matrix 
elements known is elementary and so the kinetic contribution r) k to the shear viscosity can be 
easily calculated from Eq. (|41j) . Finally, use of Eq. (|43jl in Eq. (|4*2|) determines the collisional 
transfer contribution to the shear viscosity. The result is (see Appendix EJ) 



, 2 2 

47T \ ^ x ^ 

15 

i=i j=i 



rriiTj + rrijTi 



2nm;m 



1/2 



Ci 



rrijTi 



rriiTj + rrijTi 



(50) 
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Equations (|4ip. (|46|. and (|50|) provide the explicit expression for the shear viscosity 77 of 
a dense granular binary mixture under driven USF in the first Sonine approximation. This 
coefficient is given in terms of the restitution coefficients «u, a 22 , and a± 2 , the temperature 
T, and the parameters of the mixture, namely, the masses m;, the sizes <7j, the mole fractions 
%i and the solid volume fraction <fi = <p 1 + 2 . Here, 4>i — ( 7l /^) n i a i is the species volume 
fraction of the component %. To get the explicit dependence of 77 on 0, the form of the 
pair correlation function Xij at contact must be chosen. A good approximation for Xij f° r a 



mixture of hard spheres is given by the generalized Carnahan-Starling form 



36] 



Xl] l-0 + 2(l-0)2 ^ + 2(1-0)3^J ' ^ 

where /3 = ir^nicrf + ri2<7f)/6. 

Before studying the general dependence of 77 on the parameter space, let us consider some 
special limit cases. In the elastic limit, an = a 22 = c*i2 = 1, C — 0> li — 1> $i — l/V^i, 
# 2 = I///12, and ci = c 2 = 0. In this case, the shear viscosity coefficient of an unforced 
(£ = 0) mixture can be written as 



V 



El 877^^ , \ k 4 v-^ / 1 r nmimjT\ 1 ^' 2 , 

1 + 7F a >j\>j"Jl l J' U + ( m . + m . J ( 52 ) 

i=l V j=l / i=l i=l ^ J ' 



where now the kinetic contributions r/f verify the set of equations (jUjj) with = = 
and 

16 A nec&Xit ( 2ifTm t \ 1/2 
^•=15^ (^ + ^)3/2 {-^-) (SrnA j +3m i 5 ij -2m i 6 j t), (53) 



T 8^71^- 

Ay = -^-zr-HjiXij- (54) 



^ 3 

15 T 

Equations (|52j) - (|54j) agree with the first Sonine a ppr oximation to the coefficient of shear vis- 
cosity of a molecular gas-mixture of hard spheres |37|]. In the case of mechanically equivalent 
(inelastic) particles, 7i = 1, Ci = C2 = C an d c\ = c 2 = c, where 
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FIG. 1: Plot of the reduced shear viscosity rf as a function of the restitution coefficient a for a 
binary mixture with parameters <j) = 0, x\ = 1/2, (X\jai = 1, and m\jm2 = 4 in (a) the unforced 
case (£ (0) = 0) and (b) the forced case = C {0) )- 



In this case, Eqs. flUJ), (gUJ), and fl5UJ) yield 



7] = Tj" 1 + 

and the kinetic part r? k is 



40 X (l + «) 



nT 



V 



where 



^-^ (0) + C (0) ) 

16 o /ttT 



4 /nT c 
+ 5 Y — <MX^(1 + «)(1 - gj), 



l--(l + «)(l-3a)0 X 



-ncr 



-X 



1 - -(1 

4 V 
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(57) 



(58) 



(59) 



5 \ m 

The expression (J57|) coincides with the one recently obtained for a granular monocomponent 
gas 0,1^. Finally, when — > it is easy to check that the results derived here reduce to 
those previously found in Ref. for a dilute gas. This shows the self-consistency of the 
present description. 

Before comparing the kinetic theory predictions with numerical simulation data, it is 
instructive to compare the results obtained in the unforced (^°- ) = 0) and driven = 
cases. In Fig. ^ we plot the reduced shear viscosity rj* as a function of the (common) 
restitution coefficient = a for 01/02 = 1, mxjm^ = 4, X\ = 1/2, and = in the above 
two cases. Here, the reduced shear viscosity rf is defined as 

v 



n 



nT 



(60) 
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where 

v = y/nnaf 2 v (61) 

is an effective collision frequency. We see that the Navier-Stokes shear viscosity of the 
(unforced) gas differs from the shear viscosity of the gas when the latter is excited by the 
(Gaussian) external force, the discrepancy increasing as the restitution coefficient decreases. 
This shows again that the driving force does not play a neutral role in the problem and 
the transport property is affected by this type of external forcing mechanism [6]. However, 
for practical purposes, the introduction of these driving forces has the advantage of that 
they can be incorporated into the kinetic theory very easily and they allow, for instance, to 
test the validity of some of the underlying assumptions made in the theory through a direct 
comparison with computer simulations. 



V. MONTE CARLO SIMULATION FOR UNIFORM SHEAR FLOW 

The expression obtained in the previous section for the shear viscosity requires the trun- 
cation of an expansion of the integral equations in Sonine polynomials. To assess the degree 
of accuracy of this approximation, one has to resort to numerical solutions of the Enskog 
equation, such as those obtained from Monte Carlo simulations. In this section, we briefly 
describe the method employed in the simulation in the case of the USF state. 

For a granular fluid under USF and in the absence of a thermostating force (£ = 0), 
the energy balance (fTHjl leads to a steady state when the viscous heating effect is exactly 
balanced by the collisional cooling[38]. However, when the granular mixture is excited by 
the Gaussian force 

Ft = ~m^V, (62) 

that exactly compensates for the collisional energy loss (£ = £), the viscous heating domi- 
nates and the temperature obeys the equation 

dT 2 

-dt = ~Yn aP ^ (63) 
Since the granular temperature T increases in time, so does the collision frequency u(t) oc 
y/T(t), and hence the reduced shear rate a*(t) = a/u(t) (which is the relevant uniformity 
parameter) monotonically decreases in time. Under these conditions, the system asymptot- 
ically reaches a regime described by linear hydrodynamics and the (reduced) Navier-Stokes 
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shear viscosity rf can be measured as 

?f = -lim% (64) 

where P* y = P xy /nT. Recently, this idea has been used to identify the shear viscosity of 
a (heated) granular binary mixture in the low-density regime [3]. The comparison with 
kinetic theory showed an excellent agreement over a wide range of values of the restitution 
coefficient and the rest of parameters characterizing the system. 

We have numerically solved Eq. (J2(Jj) by means of an extension of the well-known Direct 
Simulation Monte Carlo (DSMC) method [19] to dense gases. The method is usually referred 
to as the Enskog Simulation Monte Carlo (ESMC) method. This method was devised to 
mimic the dynamics involved in the Enskog collision term and it has been previously used 



to analyze rheological properties of a elastic dense gas 23] and the shock-wave structure 39]. 
In the present work, the ESMC algorithm has been modified to study the dynamics of a 
granular binary mixture of a finite density. Since the USF is spatially homogeneous in the 
local Lagrangian frame, the simulation method becomes especially easy to implement and 
efficient from a computational point of view. This is an important advantage with respect 
to molecular dynamics simulations. Nevertheless, the restriction to this homogeneous state 
prevents us from analyzing the possible instability of USF or the formation of clusters or 
microstructures. 

The ESMC method as applied to a granular binary mixture under USF is as follows. The 
velocity distribution function of the species i is represented by the peculiar velocities {Vfc} 
of Ni "simulated" particles: 

, ^ 

fi (y,t) ^nt-^SiV -V k (t)) . (65) 
1 k=i 

Note that the number of particles Ni is arbitrary, but must be taken according to the relation 
N\/N 2 = ni/n 2 - At the initial state, one assigns velocities to the particles drawn from the 
Maxwell-Boltzmann probability distribution: 

fity, 0) = n % vr" 3 U T 3 (0) exp (-V 2 /V ](0)) , (66) 

where V^(0) = 2T(0)/mj and T(0) is the initial temperature. To enforce a vanishing initial 
total momentum, the velocity of every particle is subsequently subtracted by the amount 
A^ -1 J2k Vfc(0). In the ESMC method, the free motion and the collisions are uncoupled 
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over a time step At which is small compared with both the mean free time and the inverse 
shear rate. As a* decreases monotonically in time, the value of At must be updated in the 
course of the simulation. In the local Lagrangian frame, particles of each species (i = 1, 2) 
are subjected to the action of a non- conservative inertial force F$ = — rrii a • V. Thus, 
the free motion stage consists of making — > — a • V^At. In the collision stage, 
binary interactions between particles of species % and j must be considered. To simulate 
the collisions between particles of species i with j a sample of ^NiUmaxAt pairs is chosen 
at random with equiprobability. Here, Wmax is an upper bound estimate of the probability 
that a particle of the species i collides with a particle of the species j. Let us consider a 
pair (k, £) belonging to this sample. Hereafter, k denotes a particle of species i and i a 
particle of species j. For each pair (k,£) with velocities (V&, Vg), the following steps are 
taken: (1) a given direction &ki is chosen at random with equiprobability; (2) the collision 
between particles k and t is accepted with a probability equal to &(gke • ^ke)^kf V^max, 
where u^f = AnafjUjlgke • &u\ an d gu — — *Vt — er^-a ■ tr^; (3) if the collision is accepted, 
postcollisional velocities are assigned to both particles according to the scattering rules: 

V fc -> V fe - (1^(1 + aiij)(gk£ ■ <Tki)Vki , (67) 



V £ -> V £ + J4y(l + Oy)(gfcl • °kl)Vkl • (68) 

If in a collision > cjmax, the estimate of cjmal is updated as Wmax = ^4^- The procedure 
described above is performed for i = 1, 2 and j = 1,2. The granular temperature is calculated 
before and after the collision stage, and thus the instantaneous value of the cooling rate ( is 
obtained. After the collisions have been calculated, the thermostat force (J52*j) is considered 
by making V fc -> V fc + 1/2 CV fe Ai 

In the course of the simulations, one evaluates the kinetic and collisional transfer contri- 
butions to the pressure tensor. They are given as 

2 Ni 

1=1 1 k=l 

y^J 1 :! 1 ".^:^ 1 + a ij)(gH • Vkt)VktVkt , (70) 



2iVAt 

ke 



where the dagger means that the summation is restricted to the accepted collisions. The 
shear viscosity r\ is obtained from (jSljl . To improve the statistics, the results are averaged 
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FIG. 2: Plot of the ratio rj*(a,<p)/r]*(l,<p) for a mono component gas as a function of the solid 
fraction 4> for two different values of the restitution coefficient a: (a) a = 0.9 (circles), and (b) 
a = 0.8 (squares). The lines are the theoretical predictions and the symbols refer to the results 
obtained from Monte Carlo simulations. 

over a number M of independent realizations or replicas. In our simulations we have typically 
taken a total number of particles N = N\ + N2 = 10 5 , a number of replicas M = 10, and 
a time step At = 3 x 1CT 3 Aii/Voi(0), where An = (v^r^iOn) -1 is the mean free path for 
collisions 1-1. 

VI. RESULTS 

In this section we compare the results obtained from the Chapman- Enskog method for 
the shear viscosity coefficient of a heated granular mixture (i.e., with = with 
those obtained from the ESMC method. For the sake of simplicity, we assume that 
a n — OL22 = CC12 = a, so that we reduce the parameter set of the problem to five quan- 
tities: {a, (f), m\jm2i o\j(Ji, x{\. For concreteness, henceforth we will assume that m\ > m 2 
and 0i > (72- To compare and contrast the results of a binary mixture with that of its 
monocomponent counterpart, we first show some results for a monodisperse system over a 
range of solid fractions and restitution coefficients. 
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FIG. 3: Plot of the reduced shear viscosity rf of a mono component gas as a function of the 
restitution coefficient a for three different values of the solid fraction (ft: (a) (ft = (circles), (b) 
(ft = 0.2 (squares), and (c) (ft = 0.4 (triangles). The lines are the theoretical predictions and the 
symbols refer to the results obtained from Monte Carlo simulations. 

A. Monocomponent dense gas 

Figure El shows the dependence of the ratio r]*(a,(ft)/ri*(l,(ft) on the solid fraction (ft for 
two values of the restitution coefficient. The symbols represent the simulation data while 
the lines refer to the theoretical results obtained from the Enskog equation, Eqs. (J57j) -(|59| ) . 
Both theory and simulation show that, for a given value of the density, the shear viscosity 
increases with decreasing a (i.e., greater dissipation) if the solid fraction is smaller than a 
threshold value (fto(a), while the opposite happens if (ft > (fto(a). Similar threshold values 
exist for the kinetic and collisional parts of the shear viscosity. We observe that in the range 
0.8 < a < 1, the kinetic theory calculations show that these threshold values are practically 
independent of the restitution coefficient. Specifically, (fto(a) ~ 0.16, while the corresponding 
values for the kinetic and collisional parts are, respectively, 0.23 and 0.05. It is apparent 
that the comparison between Monte Carlo simulation data and theoretical results shows an 
excellent agreement over the entire range of densities considered. The dependence of rj* (a, (ft) 
on dissipation is plotted in Fig. El for three different values of the solid fraction. We see that 
in general the influence of dissipation on the shear viscosity i]*(a) is quite significant, except 
for (ft = 0.2 which is very close to the threshold value (ft . As in Fig. |2l the theory compares 
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FIG. 4: Plot of the kinetic part 77 of the reduced shear viscosity as a function of the solid fraction 
<j) for m,\jmi = 4, o\joi = 1, x\ = 1/2 and three different values of the restitution coefficient a: (a) 
a = 0.9 (circles), (b) a = 0.8 (squares), and (c) a = 0.7 (triangles). The lines are the theoretical 
predictions and the symbols refer to the results obtained from Monte Carlo simulations. 

quite well with simulation data, except perhaps at <fi = 0.4 for strong dissipation (a = 0.6). 
B. Binary dense mixture 

Now, we consider granular binary mixtures whose particles can differ in size and mass. 
First, to analyze density effects on the shear viscosity, in Figs. 0] and 03 the parameters 
of the mixture are mi/m 2 = 4, o\j<J2 = 1, and x\ = 1/2. Three different values of a 
are studied: a = 0.9, 0.8, and 0.7. The symbols are the same as in the previous figures. 
Figure 0] shows the dependence of the kinetic part r] k * = vr^jnT on the solid fraction 0, 
while the total shear viscosity rj* is plotted in Fig. The good agreement between theory 
and simulation indicates that both kinetic and collisional transfer contributions are given 
accurately by the first Sonine approximation. As in the monocomponent case [cf . Fig. |2] , the 
shear viscosity of a granular mixture decreases (increases) as the inelasticity increases if the 
solid fraction is larger (smaller) than a given threshold value <po- The value of 0o depends 
on the parameters of the mixture although it is practically independent of dissipation. For 
the mixture considered in Fig. El (po(a) ~ 0.22. 

Next, we explore the influence of dissipation on the reduced shear viscosity 77* for different 
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FIG. 5: Plot of the reduced shear viscosity rj* as a function of the solid fraction 4> for m\/m2 = 4, 
C1/C2 = I, and x\ = 1/2 and three different values of the restitution coefficient a: a = 0.9 (solid 
line and circles), a = 0.8 (dashed line and squares), and a = 0.7 (dotted line and triangles). The 
lines are the theoretical predictions and the symbols refer to the results obtained from Monte Carlo 
simulations. 
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FIG. 6: Plot of the reduced shear viscosity r]* as a function of the mass ratio m\/m2 =, for <\> = 0.2, 
ci/c2 = 1, xi = 1/2 and three different values of the restitution coefficient a: a = 0.9 (solid line 
and circles), a = 0.8 (dashed line and squares), and a = 0.7 (dotted line and triangles). The lines 
are the theoretical predictions and the symbols refer to the results obtained from Monte Carlo 
simulations. 
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FIG. 7: Plot of the reduced shear viscosity rj* as a function of the size ratio <xi/o"2 for (f) = 0.2, 
mi/rri2 = 4, x\ = 1/2 and two different values of the restitution coefficient a: a = 0.9 (solid line 
and circles) and a = 0.7 (dashed line and triangles). The lines are the theoretical predictions and 
the symbols refer to the results obtained from Monte Carlo simulations. 

values of the mass ratio, the size ratio, and the mole fraction. We consider a solid fraction 
<p = 0.2 and different values of the restitution coefficient. In Fig. El we plot rj* versus the 



mass ratio m\jmi for <J\/<Ji = 1 and x\ — 1/2. As in the low-density case 181 ] . we see that 
the influence of dissipation on rj* becomes important as the mass disparity increases. At a 
given value of the mass ratio, rf decreases (increases) with dissipation if the mass ratio is 
smaller (larger) than a certain threshold value, which value seems to be again practically 
independent of the restitution coefficient. Regarding the comparison between kinetic theory 
and simulation, we see that the agreement between both approaches is similar to the one 
previously obtained, although the discrepancies tend to increase as a decreases. Figure [7| 
shows the results for rf as a function of the size ratio for m\/rri2 = 4 and x\ = 1/2. We 
observe that the influence of a on 77* is less significant as the one found before in Fig. El for 
the mass ratio. Finally, in Fig. [HI V* is plotted as a function of the concentration ratio X\jxi 
for 777.1/777.2 = 4 and a\ja% = \. As in Fig. [7J theory and simulation predict a weak influence 
of dissipation on the shear viscosity over the range of values of composition considered. It 
is worthwhile noting that the trends observed in Fig s. EHU for finite density are similar to 
those previously reported in the low-density limit [18]. 
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FIG. 8: Plot of the reduced shear viscosity r/* as a function of the concentration ratio xi/x% for 
<j> = 0.2, m\/m2 = 4, o\jo2 = 1 and two different values of the restitution coefficient a: a = 0.9 
(solid line and circles) and a = 0.7 (dashed line and triangles). The lines are the theoretical 
predictions and the symbols refer to the results obtained from Monte Carlo simulations. 

VII. DISCUSSION 

The main goal of this paper has been to determine the shear viscosity rj of a binary 
mixture of smooth inelastic hard spheres described by the Enskog equation. To get the 
dependence of r\ on the parameters of the mixture, the special state of uniform shear flow 
(USF) has been considered. The USF is characterized by constant partial densities, uniform 
temperature, and by a linear profile of the x-component of the flow velocity along the 
y-direction. The (constant) shear rate a is the relevant nonequilibrium parameter of the 
problem. Two complementary approaches have been used. First, a normal solution to the 
Enskog equation is obtained through first order in a by means of the Chapman-Enskog 
method. As in the elastic case, the shear viscosity coefficient rj is given in terms of the 
solution of a set of coupled linear integral equations, which are solved approximately by 
taking the leading terms in a Sonine polynomial expansion. The explicit form of r\ is given 
by Eqs. pH ]) -(|3U jl as a function of the restitution coefficients, the temperature, the total 
solid fraction, and the masses, sizes, and concentrations of the constituents of the granular 
mixture. Second, the Enskog equation has been numerically solved in the USF by using 
an extension of the well-known Direct Simulation Monte Carlo Method (DSMC) ^| of 
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the Boltzmann equation. The simulation has been performed by introducing an external 
(Gaussian) force which heats the system to compensate for the energy lost in collisions. Due 
to the action of this external driving force, the shearing work still heats the mixture and 
so the reduced shear rate a*(t) goes to zero for long times. As a consequence, the system 
reaches a regime described by linear hydrodynamics and the shear-viscosity coefficient can 
be measured in the simulations. 

The analysis made here extends previous results ^| obtained by the authors for a binary 
mixture at low-density. As in the latter case, the comparison between the Chapman- Enskog 
results in the first Sonine approximation and simulation data shows in general an excellent 
agreement for a wide range of values of densities, dissipation and parameters of the mixture. 
Discrepancies with simulation results are due mainly to the approximations carried out in the 
Chapman-Enskog scheme, and more specifically in taking only the first Sonine correction. 
However, apart from this source of slight discrepancy, the good agreement obtained here is a 
further testimony to the validity of a hydrodynamic description for granular media beyond 
the weak dissipation limit. Moreover, a test of the utility of the Enskog theory at high 
densities is possible using molecular dynamics simulations. Previous comparisons at the 
level of partial temperatures and self-diffusion coefficient j3| indicate that the range of 
densities for which the RET applies decreases with increasing dissipation. We hope that 
the present results stimulate the performance of such simulations in the case of the shear 
viscosity. 

As in the low-density case, theory and simulation show that the dependence of 77 on 
dissipation increases as the mass differences increase. The dependence of the shear viscosity 
on inelasticity is not significantly affected when composition and diameters are changed. 
With respect to the dependence on density, the results indicate that the shear viscosity 
of the granular fluid is larger than the one corresponding to a molecular fluid if the solid 
fraction is smaller than a threshold value </>o, while the opposite happens when <fi > 0o- 
The value of 0o depends on the mechanical parameters of the mixture, but is practically 
independent of the restitution coefficients. 

Recently, a seemingly similar analysis on rheology of bidisperse granular mixtures has 
been carried out via event-driven simulations [40]. However, this study is addressed to the 
steady sheared state achieved when viscous heating and collisional cooling exactly cancel 
each other. Under these conditions, due to the coupling between dissipation and the shear 
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rate, the granular fluid is far away from the Navier-Stokes regime (non-Newtonian fluid), 
except when a — *■ 1. This precludes the possibility of making a comparison between the 



molecular dynamics results found in Ref. 40( for the shear viscosity with the predictions of 
the Enskog equation. 

The main limitation of the results derived here is its restriction to the uniform shear flow 
state. The extension of this study to more general hydrodynamic states for a dense binary 
mixture, e.g. those with gradients of concentrations and temperature as well, is somewhat 
complex due to the technical difficulties of the Chapman-Enskog method associated with 
the spatial dependence of the pair correlation functions considered in the RET. We plan to 
extend the results derived years ago from the RET by Lopez de Haro et al. [32] for a mixture of 
smooth elastic particles to the case of inelastic collisions. Once the complete hydrodynamic 
equations of the mixture is known, some insight could be gained into the understanding of 
phenomena very often observed in nature and experiments, such as separation or segregation. 
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APPENDIX A: COLLISIONAL TRANSFER CONTRIBUTIONS 

In this Appendix some details of the derivation of the collisional transfer contributions 
to the heat and momentum fluxes are given. First, we consider the collisional momentum 
transfer: 

2 2 „ 



I P = / ^VxmiViJg [r, v^fufj 

i=l j=l J 

= ^ ^ Qjj / dv 1 / dw 2 I daQ(a ■ g)(<r • g)mjVi 

1=1 1 = 1 J J J 



3 

X 



[<*i3 2 Xij (r, r - (Tij)fi(r, vi; t)fj(r - o ih v' 2 ; t) 
~Xij(r, r + o- -)/<(r, v x ; t)fj(r + (T ij} v 2 ; t)\ . (Al) 

Now, we change variables to integrate over and v' 2 instead of vi and v 2 in the first term of 
the right-hand side of (jAl|) . The Jacobian of the transformation is a io - and <x-g = — ciij (ct- g') . 
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Also, Vi(v^, v 2 ) = v" = vi — (J,ji(l + aiij)(T(er ■ g) and in addition, we make the change 
<r — > — a. Thus, the integral ()A1|) becomes 

2 2 



a ?> / dvi / rfv2 / d ° ^ ' " gW v i - v i) 
i=i j=i J J J 

*Xij( r , r + Vij)fi( r , vi; t)/j-(r + try, v 2 ; t). (A2) 
Since m^v" — Vi) = rrijiy^ — v 2 '), Eq. (|A2J) can be rewritten as 

I P = / rfvi / rfv2 / dvQft ■ g)(<r • g)m j (v 2 - v 2 ) 

i=i j=i J J J 

><Xii(r, r + (Jij)fi(r, v x ; ^(r + try, v 2 ; t) 

2 2 /i /i a 

= / rfvi / dv2 / ^©(^-g)(^-g)wi(vi-Vi) 

1=1 j=l J J J 

x Xij(* ~ <Tij,r)fj(r, v a ;t)/i(r - try, Vf, t), (A3) 

where in the last step we have exchanged the roles of the species i and j , which implies the 
changes <x — > — <x and Vi <-> v 2 . Combination of Eqs. (|A2|) and (|A3|) yields 

2 2 2 

Ip = / rfVl / rfV2 / <J?e ( 5 '8)( ?, 8Wv?-Vi) 

j=i i=i J J J 

x [Xij( r , r + o"ij)/»(r, vi; i)/,-(r + try, v 2 ; i) 

-Xy(r - try, r)/ 4 (r - try, vi; *)/j(r, v 2 ; t)] . (A4) 
Let Fy(ri,r 2 ) be the function 

F y (ri, r 2 ) = x« (ri, r 2 )/ i (ri)/ i (r 2 ). (A5) 

Next, we use the identity 

f 1 d 

Fij(r, r + try) - Fy(r - cry, r) = J d\— Fy (r - (1 - A) try, r + Aery) 



CTy-v/ d\F i5 {v- (1- A)try,r+Atry). (A6) 
Jo 
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Therefore, from Eqs. (jA5|) and (|A6|) . Eq. ()A4|) can be rewriten as a divergence: 
2 2 a 



i P = -v-EE^ / rfv i / dv 2 / da e(* 
i=i j=i j j j 

xmi(vi - v") / dA%ij[r - (1 - A) try, r + Aery] 
x/j(r - (1 - A)<Tij, vi;t)/j(r + Aery, v 2 ;t) 

- -V • £ £ 4^;^ / -v. / «v 2 / «B 0(9 • g)(9 ■ g) W 



i=l j = l 

x 



/ d\xij[r - (1 - A)cry,r + Acry]/^ - (1 - \)cr ij ,v 1 ;t)f j (r + Aery, v 2 ;t). 
Jo 



(A7) 

According to the momentum balance equation ©, the divergence of the collisional transfer 
part P c is 

2 2 „ 

E E / dvm ^ v " u ) J SM/<. /;] = - v ■ pc - ( A8 ) 

i=l j=l J 

By taking into account (jA8|) . one directly gets the expression (JT3j) for the collisional part of 
the pressure tensor: 

2 2 

Jmi + rrij 2 J J J 

1=1 J = l J 

x / d\xij[r - (1 - A)cry,r + A<7y]/i(r - (1 - A)cry, vi; £)/j(r + Aery, v 2 ;t). 
./o 

(A9) 

The collisional transfer contribution to the energy balance equation is 

Ie = EE / ^V^M^,/,]. (A10) 
1=1 , = 1 ^ 
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By following similar mathematical steps as in the case of momentum, one gets 

2 2 2 

^ = EEy [dvi f ^e(£- g )(£. g )^(v;' 2 - Vl ) 
i=i j=i j j j 

x [Xij( r , r + <Tij)fi( r , vi; *)/j-(r + try, v 2 ; t) 



-Xy( r , r - o-ij)/i( r - fy, vi; *)//(r, v 2 ; t)] 

2 2 

I I I \ ~ I / / / 



I 3=1 ' ' 

><XyO, r + tTij)fi(T, vi; t)£ (r + try, v 2 ; t) 



x y( v ! - v ?) J d\xij[r - (1 - A) try, r + Aoy) 
xfi(r - (1 - A)try,v 1 ;£)/ i (r + Atry, v 2 ;*) 

1=1 j=l J 

*Xij (r, r + <r y )/i(r, v x ; £)/,(r + try, v 2 ; t). (All) 



Here, we have used the identity (|A6|) and the scattering law 



m,(V " v ) = m,(v - V) - — ^-(1 - • g )'. (A12) 



- mim ^_(i- a ?.)(^ -g) 2 

The balance energy equation (JZJ) yields 

E E / rfv T (V " ^ = - V • q C - P C : Vu - ^TC, 
i=i j=i J 



(A13) 



where q° is the collisional contribution to the heat flux and £ is the cooling rate. To identify 
such quantities, we use the relation 



Y (v? - v?) = ^^(l - a%){* • g) 2 + m iNi {l + ay)(& • g)[(<x • Gy) + (3= • u)], (A14) 

with Gy = /ijjVi + /ijjV 2 , V = v — u being the peculiar velocity. Comparing Eqs. (|A11J) and 
flA13|) and taking into account Eqs. (jA7|) and (jA14j) . one can finally obtain the expressions 
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()14j) and ()15|) for q c and (. They are given by 

2 2 



q c M) = EE-I 



TOjTOj 1 + a. 



i=i i=i 



mj + to.,- 2 



/ dvi y cfv 2 y dd-Q(S ■ g)(<7 ■ g) 2 <x 



x 



(<r • Gy) + — A*ji(l - 



x / d\Xij[r - (1 - XjtTij, r + Aer^ 



x/j(r - (1 - \)(Tij, vi; t)fj(r + Xcr^, v 2 ; £) 



, y mj + to,- 2 



X 



(ct • G y ) + - (fiji - Hij) (1 - • g) 



/ d\xij[r - (1 - A)<r# , r + Xa^f^r - (1 - A)<7y, vi; *)/,-(r + Acr^-, v 2 ; t), 
Jo 

(A15) 



2 2 



CM) = ^EE^^ 1 -* / dvi / rfV2 / ^©(^-g)(^-6) 3 

i=l j=l 1 

xXii(r,r + o-i i )/ i (r,v 1 ;t)/ i (r + <r ii ,v 2 ;t). (A16) 



The second equality in Eq. (|A15|) has been obtained by exchanging the roles of species % and 
j. In the case of mechanically equivalent particles, Eqs. (|A9I) . ()A15|) . and ()A16|) reduce to 
those previously obtained for a monocomponent dense gas [21 1. 



APPENDIX B: FIRST ORDER SOLUTION TO THE USF 

In this Appendix we apply the Chapman- Enskog method to solve Eq. (|2U|) to first order 
in the shear rate a. First, in order to get the kinetic equation for the Enskog collision 
operator (|2*Tj) must be expanded as 

jf, - ^[^./fi+^'i/f'./ri+^'i/f'./] 1 '] 

+aA«[/f ) ,/f], (Bl) 
where the (Boltzmann) collision operator J^[Xi,X 2 ] is defined in Eq. ()31|) and 



AiifVrljfjf] = xaoij dv 2 J dae(a. g )(a. s )a y 



+/rcvx)^/r(v a 



<9 



(0), 



(B2) 
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Therefore, the distribution verifies the equation 



a, 



(0) 



1,(0) A 

2^ <9V 



V + A) if + Mjf = aVy ^ + a ± A, [if , /f ] 

' X 1=1 



(B3) 



where it is understood that i ^ j on the left hand side and the linear operators C{ and Aii 



arc 



A/* 



(i) 



(1) r(0)l 



r(l) /(0)- 



(B4) 



(0) ,(l)i 



(B5) 



In these equations, use has been made of the fact that cf jf = according to the second 
identity of Eq. ([28)1 . Furthermore, = by symmetry because V ■ u = in the USF. The 
action of the time derivative cf on /f can be easily obtained from (J28)) as 

9 (o) / (i) = _ (c (o)_ e (o) )T5Ti .(i) ; (B6) 

and so, the integral equation (jB3|) finally becomes 



_ {C (°)_£(o) )TdT + ^(o)A_.y + Ci 



This is the result fl38j) used in the text. 



(0) -(0)1 



(B7) 



APPENDIX C: EVALUATION OF SOME COLLISION INTEGRALS 



In this Appendix some collision integrals appearing along the text are evaluated. First, 
let us consider the integral (|49|) 



A; 



rrij 



dV l m l V lx V ly K lJ [V 1 \fl°\f. 



ill) ^(0)j 



/ rfVi / ^V 2 / d*Q{v-z){*-z)a y V lx V, 



+// 0) (Vi)^r- /f(V 



<9 



(o), 



2x 



(CI) 



As done in Appendix [XJ we change variables to integrate over and V2 instead of Vi and 
V2 in the first term of the right-hand side of (jClj) . Thus, the integral becomes 



A; 



rrii 



1 i 



d 



■2x 



(C2) 
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where 



V" = Vi - ^(1 + ay)£(£ ■ g). (C3) 
Using (jC3j) . the last term on the right hand side of (jC2|) can be explicitly computed as 



iy 



-fJ-jiO- + Oiij)(a- ■ g) [Gi j)X a y + G ijt y(j x 

- Hji{9x<? y + g y Vx) - + a>ij)(o- ■ g)a x a y ] , 



(C4) 



where Gy ;X = [iijV\ x + fijiV 2x and Gij )V = ^%jV\ y + l^j^Yiy Substitution of Eq. (|C4|) into Eq. 
2J) allows the angular integral to be performed with the result 



/2ir 
da Q(a ■ g)(5 ■ g)a y {v;' x V(' y - V lx V ly ) = - ^(1 + ay) [(2^ + g 2 )G lh2 



+2g x g y G ijt y + - 3^)^^ + -/i,i(4 - 3aij)g x g 2 



(C5) 



Using (|C5j) . the integral (|C2|) becomes 
27r m,- 



A 



15 n^ 2 



Xy4^(l + «y) / dVi I dV 2 f \v 1 )fP(V 2 ) 



f(0). 



X 



3 

27r min,; 



4 



^(3c*y - 1)(^ 2 + V 2 2 ) - - (mjV? - Hi vT) 



"%>«] 



15 T 



Ot; 



Ti T 7 - 
— + — 

m,- m, 



T- — T 



(G6) 



In the case of mechanically equivalent particles, Eq. (|C6J) coincides with the one previously 
obtained in the context of determining the shear viscosity in a monocomponent granular gas 



The collision frequencies ry defined by the integrals (}4Tj) and (|4~81) are the same as those 
appearing in the Boltzmann limit (except for the factors Xy)|l2l The details will not 
be repeated here and only the results are displayed. They are given by 



16 irTt 2 

Til = ~=-\ WiO-iXll 

5 V mi 



1-4(1 -a n ) 



2L 
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+—^2^2/121X12^0(1 + ' [60^(^ 12 9 2 - // 2 i0i)(0i + ^ 2 )" 1/2 

+^/i2i^ 2 (^i + ^ 2 ) 1/2 (3 - a l2 ) + 50 1 ~ 1 (0i + ^ 2 )- 1 / 2 



c 2 2^(12/121 + 9/ii 2 - 10) - X (5 - 6/i 21 ) - |/x 21 (3 - a 12 )(0i + 9 2 ) 



16 



(6 1 + 6 2 )^ 



(C7) 
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ri2 = ^VWj^— xWi + ^2)o'i /2 6 2 1/2 [Qe^(fi 12 d 2 - f i 21 e 1 )(e 1 + 6 2 ) 

15 /ii2 ~ 



-1/2 



+ 2^i^ 2 (^i + ^) 1/2 (3-«i 2 ; 



50," 1 



2 I 17 ! 



1-1/2 



d 20 1 (1O - I2fi 12 - 9fi 21 ) + 2 (5 - 6^ 12 ) - |// 2 i(3 - ai 2 )(0i + 2 ) 



(Ci 



16 (^i + 2 ) 5 / 2 

The corresponding expressions for t 22 and r 21 can be inferred from Eqs. (jC7|) and (|C8|) by 
exchanging 1 <->• 2. 

Finally, the collision contribution to the shear viscosity is given by Eq. (|42j) . To get 
explicit results, we have to evaluate integrals of the form 



A l3 = J dV, J dV 2 fi°\Y 1 )ff\v 2 )g, 



(C9) 

by using the leading Sonine approximation (j4"3"|). Let us consider the integral Ai 2 . Substi- 
tution of Eq. ()43|) into Eq. ()C9|) and neglecting nonlinear terms in q, A 12 can be written 

as 



A 12 = n 1 n 2 7r- 3 v (e 1 e 2 f 2 {i(e 1 ,e 2 ) 



4(^+«^ + ^| M 



d 15 
'd02~ + T 



(CIO) 



where the dimensionles integral I(6i,6 2 ) is 

J(0 ls 2 ) = fdV*if dV* 2 e- e ^ 2 - e ^\ 
with V* = V/f . The integral 7(0i, 9 2 ) can be performed by the change of variables 

x = vj- v 2 *, y = e x \i + e 2 v* 2 , 

with the Jacobian (9\ + 9 2 )~ z . The integral becomes 



(Cll) 



(C12) 



(C13) 



Use of this result in (jC10|) gives 

Ala = A21 



2 /0i + 2 

— nin 2 w 

7T V 0102 



1/2 



Cl / 2 



<~2 



16\0i + 02/ 16 \0i + 02 



01 



(C14) 
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The corresponding expressions for An and A22 can be easily inferred from Eq. (jC14j) by 
exchanging 1 <-> 2: 

4 2 [W ci 



■^^V^-sb (C15) 



^-^V^O-D- (C16) 

Equations ()C14|) - ()C16|) lead directly to the result given by Eq. (|5U|) in the text. 
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